#!/usr/bin/perl -w

# Copyright 2013 Thomas H. Schmidt
#
# This file is part of DanceSteps.
#
# DanceSteps is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# any later version.
#
# DanceSteps is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with DanceSteps; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.


### Load Packages ##############################################################
use strict;
use IO::Handle;
use FindBin qw($RealBin); # Absolute path to THIS script.
use lib ($RealBin."/modules/FileIO", $RealBin."/modules");
autoflush STDOUT 1; # For direct output.

use Commandline;
use FileIO::Gro;
use FileIO::Ndx;
use FileIO::Xtc;
use PTE qw(getElementData isElement printPte);
use MolOps qw(getComPbc);
use Grid3D;
use Statistics;
#use protein;
################################################################################



### Default Parameters #########################################################
our $version        = "rc1";              # Version number.
our $year           = "2013";             # Year (of change).

our $verbose        = 0;                  # Be loud and noisy (and global); default: silent.

my $coordInFile     = 'system.gro';       # Input coordinates file.
my $xtcFile         = 'traj.xtc';
my $ndxInFile       = '';                 # Input GROMACS index file.
my $gridOutFile     = 'grid.gro';
my $groupResname    = '^POP[ACEG]?$';       # Regex for the automatic detection of lipid molecules in the membrane GRO file.


my %timeRange;
@{$timeRange{'dump'}} = (-1);
$timeRange{'range'}{'min'}  = 0;
$timeRange{'range'}{'max'}  = 0;
$timeRange{'range'}{'step'} = 0;

my $gridDeltaX    = 0.03;                 # Grid size in x direction [nm].
my $gridDeltaY    = 0.03;                 # Grid size in y direction [nm].
my $gridDeltaZ    = 0.03;                 # Grid size in z direction [nm].
my $alongAxis     = 'z';                  # Along this axis the measurement is performed.
my $map2Slice     = 0;                    # Map all grid points to one slice and calculate the area.
my $pbcWrap       = 0;                    # Wrap the coordinates at the edges of the box.
my $vsDist        = 0;                    # Measure the thickness as a function of distance to another group of atoms.
my $radStepwidth  = 0.3;                  # Stepwidth for measuring the radial thickness (vsDist option).
my $helpAndQuit   = 0;                    # Print out program help.
################################################################################



### Internal parameters ########################################################
my %coordData;         # Filled by "GROFiles::readGro(<GROFILE>)".
my @ndxData;           # Filled by "NDXFiles::readNdx(<NDXFILE>)".
my @thickNdxGroupIds;  # Filled by "NDXFiles::selectGroupIds(...)".
my @vsDistNdxGroupIds; # Filled by "NDXFiles::selectGroupIds(...)".
################################################################################



### Print out program headlines ################################################
printHead();
################################################################################



### Handle commandline parameters ##############################################
addCmdlParam('scalar', 's',       'Input',       \$coordInFile,                $coordInFile, 'Structure file: gro');
#addCmdlParam('scalar', 'f',       'Input, Opt.', \$xtcFile,                    $xtcFile, 'Trajectory file: xtc');
addCmdlParam('scalar', 'n',       'Input, Opt.', \$ndxInFile,                  $ndxInFile, 'Index file');
addCmdlParam('scalar', 'o',       'Output',      \$gridOutFile,                $gridOutFile, 'Structure file: gro');
addCmdlParam('flag',   'h',       'bool',        \$helpAndQuit,                $helpAndQuit ? 'yes' : 'no', 'Print help info and quit');
addCmdlParam('scalar', 'b',       'time',        \$timeRange{'range'}{'min'},  ($timeRange{'range'}{'min'} ? $timeRange{'range'}{'min'} : 'no'), 'First frame (ps) to read from trajectory');
addCmdlParam('scalar', 'e',       'time',        \$timeRange{'range'}{'max'},  $timeRange{'range'}{'max'} ? $timeRange{'range'}{'max'} : 0, 'Last frame (ps) to read from trajectory');
addCmdlParam('scalar', 'dt',      'time',        \$timeRange{'range'}{'step'}, $timeRange{'range'}{'step'} ? $timeRange{'range'}{'step'} : 0, '');
addCmdlParam('array',  'dump',    'time',        \@{$timeRange{'dump'}},       $timeRange{'dump'} ? join(" ", @{$timeRange{'dump'}}) : '-1', 'Dump frame nearest specified time (ps)');
addCmdlParam('scalar', 'dx',      'real',        \$gridDeltaX,                 $gridDeltaX, 'Grid spacing along the X axis');
addCmdlParam('scalar', 'dy',      'real',        \$gridDeltaY,                 $gridDeltaY, 'Grid spacing along the Y axis');
addCmdlParam('scalar', 'dz',      'real',        \$gridDeltaZ,                 $gridDeltaZ, 'Grid spacing along the Z axis');
addCmdlParam('scalar', 'd',       'real',        \$alongAxis,                  $alongAxis, 'Axis for distance calculation.');
addCmdlParam('flag',   'map2sl',  'real',        \$map2Slice,                  $map2Slice, 'Map all grid points to one slice and calculate the area.');
addCmdlParam('flag',   'pbcwrap', 'bool',        \$pbcWrap,                    $pbcWrap ? 'yes' : 'no', 'Wrap the coordinates at the edges of the box.');
addCmdlParam('flag',   'vsdist',  'bool',        \$vsDist,                     $vsDist ? 'yes' : 'no', 'Measure the thickness of group 1 as a function of distance to group 2 (NDX file required).');
addCmdlParam('flag',   'v',       'bool',        \$verbose,                    $verbose ? 'yes' : 'no', 'Be loud and noisy');
addCmdlParam('scalar', 'lr',      'string',      \$groupResname,               $groupResname, 'Regex for the detection of the bilayer atoms of the membrane GRO file (obsolete if a membrane NDX file is given)');

cmdlParser();
################################################################################



### Print program help if the user set the flag ################################
printHelp(getCmdlParamRef(), 1) if $helpAndQuit;
################################################################################



### Read the XTC file ##########################################################
#my $xtcDataRef = FileIO::XTC::readXtc($xtcFile, \%timeRange);
#
#
#printf("\n\nFrame Step      Time     Atom cooX     cooY     cooZ\n");
#for (my $fr=0; $fr<@{$xtcDataRef}; $fr++) {
#    for (my $i=0; $i<1; $i++) {
#       printf "  %3d % 9d % 8.3f % 4d % 8.3f % 8.3f % 8.3f\n", $fr, $$xtcDataRef[$fr]{'step'}, $$xtcDataRef[$fr]{'time'}, $i + 1, $$xtcDataRef[$fr]{'atoms'}[$i]{'cooX'}, $$xtcDataRef[$fr]{'atoms'}[$i]{'cooY'}, $$xtcDataRef[$fr]{'atoms'}[$i]{'cooZ'};
#    }
#}
#exit;
################################################################################


### Read the GRO file ##########################################################
%coordData = FileIO::GRO::readGro($coordInFile); # Read input GRO file.
################################################################################


### Read the NDX file ##########################################################
if ($ndxInFile) {
    @ndxData = FileIO::NDX::readNdx($ndxInFile); # Read input NDX file.
    die "ERROR: Cannot find NDX data.\n" unless (@ndxData);
    FileIO::NDX::printNdxGroups(@ndxData);
    @thickNdxGroupIds  = FileIO::NDX::selectGroupIds(\@ndxData, 'thickness measurement', 1); # Select only one group.
    @vsDistNdxGroupIds = FileIO::NDX::selectGroupIds(\@ndxData, 'measuring the thickness within a radius to this group', 1) if $vsDist; # Select only one group.
}
elsif ($groupResname) {
    print "No NDX file defined. Will use regex of residue names to select atoms for calculation.\n";

    $thickNdxGroupIds[0] = 1;
    $ndxData[$thickNdxGroupIds[0]]{'groupName'} = $groupResname;
    for (my $i=1; $i<@{$coordData{'atoms'}}; $i++) {
        push(@{$ndxData[$thickNdxGroupIds[0]]{'atoms'}}, $i) if $coordData{'atoms'}[$i]{'resName'} =~ /$groupResname/;
#        print "Found " . scalar(@{$ndxData[$thickNdxGroupIds[0]]{'atoms'}}) . " atoms of known bilayer resnames\n";
    }
}
else {
    printHelp();
}
################################################################################



if ($vsDist) {
    print "Measure the thickness within a radius\n";

    ### Add atom masses from PTE to the coordinate data ########################
    addAtomMasses($ndxData[$vsDistNdxGroupIds[0]]{'atoms'}, $coordData{'atoms'});
    ############################################################################


    ### Detect the center of mass of group @vsDistNdxGroupIds ##################
    my %com = getComPbc($ndxData[$vsDistNdxGroupIds[0]]{'atoms'}, $coordData{'atoms'}, $coordData{'box'});
    printf("  Center of mass of group % 2d: %8.3f %8.3f %8.3f\n", $vsDistNdxGroupIds[0], $com{'cooX'}, $com{'cooY'}, $com{'cooZ'});
    ############################################################################


    ### Map group one to a grid based on the com of group 2 ####################
    # Translate each atom: a' = a + t; t = c_box - com
    # Build the grid using the existing routines
#    $radStepwidth;
    ############################################################################

    exit;
}



sub addAtomMasses {
    my $atomIdsRef = shift;
    my $coordDataRef = shift;

    foreach (@{$atomIdsRef}) {
        $$coordDataRef[$_]{'mass'} = getStdWeightByElement(atomName2Element($$coordDataRef[$_]{'atomName'}));
    }
}



sub atomName2Element {
    my $atomName = shift;

    return $atomName if isElement($atomName);

    $atomName =~ /([A-Za-z]{0,3})/;
    return $1 if isElement($1);

    my $oneLetter = substr($1, 0, 1);
    return $oneLetter if isElement($oneLetter);

    return;
}

### Map atoms to the 3D grid ###################################################
setGridDelta($gridDeltaX, $gridDeltaY, $gridDeltaZ); # Set the grid spacing in x, y & z.
atoms2Grid($ndxData[$thickNdxGroupIds[0]]{'atoms'}, \%coordData, 0, $pbcWrap);
################################################################################


my $completeGrid = 0; # Fill up all (also undefined) grid cells until the box limits are reached.


### Map atoms within a range to one slice ###################################### # From squaredance.
if ($map2Slice) {
    my $zSliceRef = mapGrid2Slice(getGridOccRef(), "z");

    my $sliceGroDataRef = slice2coordData($zSliceRef, "z", getGridZMin());
    FileIO::GRO::writeGro($gridOutFile, $sliceGroDataRef); # Print slice-coordinates on a 3D grid.

    my $sliceArea = getSliceArea($zSliceRef, $gridDeltaX, $gridDeltaY);
    printf("Area = %8.3f nm²\n", $sliceArea);
#    printf("Percent of the XY plane (A = %0.3f nm²) = %0.3f%%\n", $groData{'box'}{'cooX'}*$groData{'box'}{'cooY'}, ($nSliceVoxel*$gridDeltaX*$gridDeltaY*100)/($groData{'box'}{'cooX'}*$groData{'box'}{'cooY'}));
    exit;
}
################################################################################



### Write the grid to a GRO file ###############################################
if ($gridOutFile) {
    my $gridGroDataRef = grid2CoordData(getGridOccRef(), $completeGrid);
    FileIO::GRO::writeGro($gridOutFile, $gridGroDataRef); # Print grid-coordinates.
}
################################################################################


### Write out the normal profile of occupied grid slices #######################
my $xSecAreaOutFile = "xsecarea.dat";
my $zProfileRef = getZProfile(getGridOccRef());
xSecArea2File($zProfileRef, $xSecAreaOutFile, $gridDeltaX, $gridDeltaY, $gridDeltaZ);
################################################################################


### Write out the 2D thickness profile #########################################
my $profileOutFile = "thickness2d.dat";
my $profile2DRef = gridDistProfile2D(getGridOccRef());
gridDistProfile2D2File($profile2DRef, $profileOutFile, $gridDeltaX, $gridDeltaY, $gridDeltaZ);
################################################################################


### Print out statistics #######################################################
my @allValues = matrix2D2Vector($profile2DRef);
my @statistics = Statistics::all(\@allValues);
print $statistics[0]*$gridDeltaZ . " = mean\n";
print $statistics[1]*$gridDeltaZ*$gridDeltaZ . " = variance\n";
print $statistics[2]*$gridDeltaZ . " = stddev\n";
################################################################################


print "Calculate the volume of the group:\n";
my $volume = getGridVolume(getGridOccRef());
print "Volume = $volume nm³\n";



sub mapGrid2Slice {
    my $gridRef = shift;
    my $vector  = shift;
    $vector = ($vector eq "x" || $vector eq "y") ? $vector : "z";
    my @slice;

    if ($vector eq "z") {
        for (my $z=0; $z<@$gridRef; $z++) {
            next unless $$gridRef[$z];
            for (my $x=0; $x<@{$$gridRef[$z]}; $x++) {
                next unless $$gridRef[$z][$x];
                for (my $y=0; $y<@{$$gridRef[$z][$x]}; $y++) {
                    next unless $$gridRef[$z][$x][$y];
                    $slice[$x][$y] = $$gridRef[$z][$x][$y];
                }
            }
        }
    }
    elsif ($vector eq "x") {
        for (my $z=0; $z<@$gridRef; $z++) {
            next unless $$gridRef[$z];
            for (my $x=0; $x<@{$$gridRef[$z]}; $x++) {
                next unless $$gridRef[$z][$x];
                for (my $y=0; $y<@{$$gridRef[$z][$x]}; $y++) {
                    next unless $$gridRef[$z][$x][$y];
                    $slice[$y][$z] = $$gridRef[$z][$x][$y];
                }
            }
        }
    }
    elsif ($vector eq "y") {
        for (my $z=0; $z<@$gridRef; $z++) {
            next unless $$gridRef[$z];
            for (my $x=0; $x<@{$$gridRef[$z]}; $x++) {
                next unless $$gridRef[$z][$x];
                for (my $y=0; $y<@{$$gridRef[$z][$x]}; $y++) {
                    next unless $$gridRef[$z][$x][$y];
                    $slice[$x][$z] = $$gridRef[$z][$x][$y];
                }
            }
        }
    }
    return \@slice;
}



sub getSliceArea {
    my $sliceRef   = shift;
    my $gridDelta1 = shift;
    my $gridDelta2 = shift;
    my $nPixel     = 0;

    for (my $x=0; $x<@{$sliceRef}; $x++) {
        next unless $$sliceRef[$x];
        for (my $y=0; $y<@{$$sliceRef[$x]}; $y++) {
            next unless $$sliceRef[$x][$y];
            $nPixel++;
        }
    }

    return($nPixel * $gridDelta1 * $gridDelta2);
}



sub slice2coordData {
    my $sliceRef     = shift;
    my $vector       = shift;
    my $sliceVal     = shift;
    $vector = ($vector eq "x" || $vector eq "y") ? $vector : "z";
    my @sliceOnGrid;
    my $voxId   = 0;

    if ($vector eq "z") {
        for (my $x=0; $x<@{$sliceRef}; $x++) {
            next unless $$sliceRef[$x];
            for (my $y=0; $y<@{$$sliceRef[$x]}; $y++) {
                next unless $$sliceRef[$x][$y];
                $sliceOnGrid[$sliceVal][$x][$y] = $$sliceRef[$x][$y];
            }
        }
    }
    elsif ($vector eq "x") {
        for (my $y=0; $y<@{$sliceRef}; $y++) {
            next unless $$sliceRef[$y];
            for (my $z=0; $z<@{$$sliceRef[$y]}; $z++) {
                next unless $$sliceRef[$y][$z];
                $sliceOnGrid[$z][$sliceVal][$y] = $$sliceRef[$y][$z];
            }
        }
    }
    elsif ($vector eq "y") {
        for (my $x=0; $x<@{$sliceRef}; $x++) {
            next unless $$sliceRef[$x];
            for (my $z=0; $z<@{$$sliceRef[$x]}; $z++) {
                next unless $$sliceRef[$x][$z];
                $sliceOnGrid[$z][$x][$sliceVal] = $$sliceRef[$x][$z];
            }
        }
    }

    return grid2CoordData(\@sliceOnGrid);
}



sub matrix2D2Vector {
    my $matrix2DRef = shift;
    my @vector;

    for (my $x=0; $x<@{$matrix2DRef}; $x++) {
        next unless $$matrix2DRef[$x];
        for (my $y=0; $y<@{$$matrix2DRef[$x]}; $y++) {
            next unless $$matrix2DRef[$x][$y];
            push(@vector, $$matrix2DRef[$x][$y]);
        }
    }

    return unless scalar(@vector);
    return @vector;
}



sub gridDistProfile2D {
    my $gridRef = shift;
    my @minZ;
    my @maxZ;
    my @profile2D;

    ### Get the minimum and maximum per z row ##################################
    for (my $z=0; $z<@$gridRef; $z++) {
        next unless $$gridRef[$z];
        for (my $x=0; $x<@{$$gridRef[$z]}; $x++) {
            next unless $$gridRef[$z][$x];
            for (my $y=0; $y<@{$$gridRef[$z][$x]}; $y++) {
                next unless $$gridRef[$z][$x][$y];
                $maxZ[$x][$y] = $z;
                next if $minZ[$x][$y];
                $minZ[$x][$y] = $z;
            }
        }
    }
    ############################################################################


    for (my $x=0; $x<@minZ; $x++) {
        next unless $minZ[$x];
        for (my $y=0; $y<@{$minZ[$x]}; $y++) {
            next unless $minZ[$x][$y];
            next unless $maxZ[$x][$y];
            $profile2D[$x][$y] = $maxZ[$x][$y] - $minZ[$x][$y];
        }
    }

    return \@profile2D;
}



sub gridDistProfile2D2File {
    my $profile2DRef   = shift;
    my $profileOutFile = shift;
    my $gridDeltaX     = shift;
    my $gridDeltaY     = shift;
    my $gridDeltaZ     = shift;

    open(PROF2D, ">$profileOutFile") || die "ERROR: Cannot open 2D profile output file \"$profileOutFile\": $!\n";
    for (my $x=0; $x<@$profile2DRef; $x++) {
#        print scalar(@{$$profile2DRef[$x]}) . "\n"; exit;
        next unless $$profile2DRef[$x];
        for (my $y=0; $y<@{$$profile2DRef[$x]}; $y++) {
#            print $$profile2DRef[$x][$y] . " :: $x $y\n";
            next unless $$profile2DRef[$x][$y];
            printf PROF2D ("%.3f  %.3f  %.3f\n", $x*$gridDeltaX, $y*$gridDeltaY, $$profile2DRef[$x][$y]*$gridDeltaZ);
        }
    }
    close(PROF2D);
}



sub getZProfile {
    my $gridRef    = shift;
    my @profile;

    for (my $z=0; $z<@$gridRef; $z++) {
        next unless $$gridRef[$z];
        $profile[$z] = 0;
        for (my $x=0; $x<@{$$gridRef[$z]}; $x++) {
            next unless $$gridRef[$z][$x];
            for (my $y=0; $y<@{$$gridRef[$z][$x]}; $y++) {
                next unless $$gridRef[$z][$x][$y];
                $profile[$z]++;
            }
        }
    }

    return \@profile;
}



sub xSecArea2File {
    my $profileRef      = shift;
    my $xSecAreaOutFile = shift;
    my $gridDeltaX      = shift;
    my $gridDeltaY      = shift;
    my $gridDeltaZ      = shift;

    open(XSECPROF, ">$xSecAreaOutFile") || die "ERROR: Cannot open cross-sectional area output file \"$xSecAreaOutFile\": $!\n";
    for (my $z=0; $z<=@$profileRef; $z++) {
        next unless $$profileRef[$z];
        printf XSECPROF ("%.3f  %.3f\n", $z*$gridDeltaZ, $$profileRef[$z]*$gridDeltaX*$gridDeltaY);
    }
    close(XSECPROF);
}

#my @cavGeoCenterXy;

### Basic analysis #############################################################
#foreach (@inflGroupIds) {
#    printf("\n  ---------------------------------\n  Analyze group \"%s\" (ID=%d)...\r", $ndxData[$_]{'groupName'}, $_);
#
#    ### Basic analysis of the bilayer ##########################################
#    %{$lipidData[$_]} = Lipid::analyzeBilayer($ndxData[$_]{'atoms'}, $coordData{'atoms'}, $refAtomNames{'lipid'});
#    printf("\n    Upper leaflet: %d lipids\n    Lower leaflet: %d lipids\n",
#            $lipidData[$_]{'leaflets'}{'upper'}{'nLipids'},
#            $lipidData[$_]{'leaflets'}{'lower'}{'nLipids'}) if $verbose;
#    ############################################################################
#
#
#    ### Determine the XY area occupied by protein ##############################
#    printf("    Detecting protein area (upper leaflet)...\n") if $verbose;
#    my @tmpArray;
#
#    @tmpArray = Protein::analyzeProtein(\@ndxData,
#                                        $protGroupId,
#                                        \@suGroupIds,
#                                        $coordData{'atoms'},
#                                        $coordData{'box'},
#                                        $gridDelta,
#                                        ($lipidData[$_]{'leaflets'}{'upper'}{'headCooZAvg'} - 0.5),
#                                        ($lipidData[$_]{'leaflets'}{'upper'}{'headCooZAvg'} + 0.5),
#                                        $extCavDetection,
#                                        $ndxData[$_]{'groupName'} . ".upper.xy");
#    $gridAProt[$_]{'upper'} = $tmpArray[0];
#    $gridACav[$_]{'upper'}  = $tmpArray[1];
#    $gridAProt[$_]{'upper'} += $gridACav[$_]{'upper'} unless $lipidFilledCav; # If the detected cavity is empty (not lipid-filled), the excluded area for lipid packing is AProtein + ALipid.
#
#    printf("      Detected area excludable for lipid packing: %7.4f nm^2\n", $gridAProt[$_]{'upper'}) if $verbose;
#
#    if ($lipidFilledCav) { # Calculate the geometrical XY center of the cavity (upper leaflet).
#        my $tmpCounter = 0;
#        for (my $x=0; $x<=@{$tmpArray[2]}; $x++) {
#            next unless ${$tmpArray[2]}[$x];
#            print "$x\r";
#            for (my $y=0; $y<=@{$tmpArray[2][$x]}; $y++) {
#                next unless ${$tmpArray[2]}[$x][$y];
#                next unless ${$tmpArray[2]}[$x][$y] == 2; # Next loop if this grid point is not defined as cavity.
#                $cavGeoCenterXy[$_]{'upper'}{'cooX'} += $x * $gridDelta;
#                $cavGeoCenterXy[$_]{'upper'}{'cooY'} += $y * $gridDelta;
#                $tmpCounter++;
#            }
#        }
#        $cavGeoCenterXy[$_]{'upper'}{'cooX'} /= $tmpCounter;
#        $cavGeoCenterXy[$_]{'upper'}{'cooY'} /= $tmpCounter;
#        $cavGeoCenterXy[$_]{'upper'}{'cooZ'} = $lipidData[$_]{'leaflets'}{'upper'}{'headCooZAvg'};
#    }
#
#
#    printf("    Detecting protein area (lower leaflet)...\n") if $verbose;
#    @tmpArray = Protein::analyzeProtein(\@ndxData,
#                                        $protGroupId,
#                                        \@suGroupIds,
#                                        $coordData{'atoms'},
#                                        $coordData{'box'},
#                                        $gridDelta,
#                                        ($lipidData[$_]{'leaflets'}{'lower'}{'headCooZAvg'} - 0.5),
#                                        ($lipidData[$_]{'leaflets'}{'lower'}{'headCooZAvg'} + 0.5),
#                                        $extCavDetection,
#                                        $ndxData[$_]{'groupName'} . ".lower.xy");
#    $gridAProt[$_]{'lower'} = $tmpArray[0];
#    $gridACav[$_]{'lower'}  = $tmpArray[1];
#    $gridAProt[$_]{'lower'} += $gridACav[$_]{'lower'} unless $lipidFilledCav;
#
#    printf("      Detected area excludable for lipid packing: %7.4f nm^2\n", $gridAProt[$_]{'lower'}) if $verbose;
#
#    if ($lipidFilledCav) { # Calculate the geometrical XY center of the cavity (lower leaflet).
#        my $tmpCounter = 0;
#        for (my $x=0; $x<=@{$tmpArray[2]}; $x++) {
#            next unless ${$tmpArray[2]}[$x];
#
#            for (my $y=0; $y<=@{$tmpArray[2][$x]}; $y++) {
#                next unless ${$tmpArray[2]}[$x][$y];
#                next unless ${$tmpArray[2]}[$x][$y] == 2; # Next loop if this grid point is not defined as cavity.
#                $cavGeoCenterXy[$_]{'lower'}{'cooX'} += $x * $gridDelta;
#                $cavGeoCenterXy[$_]{'lower'}{'cooY'} += $y * $gridDelta;
#                $tmpCounter++;
#            }
#        }
#        $cavGeoCenterXy[$_]{'lower'}{'cooX'} /= $tmpCounter;
#        $cavGeoCenterXy[$_]{'lower'}{'cooY'} /= $tmpCounter;
#        $cavGeoCenterXy[$_]{'lower'}{'cooZ'} = $lipidData[$_]{'leaflets'}{'lower'}{'headCooZAvg'};
#    }
#    ############################################################################
#
#    printf("  Analyze group \"%s\" (ID=%d): Finished\n", $ndxData[$_]{'groupName'}, $_) unless $verbose;
#    print "  ---------------------------------\n\n";
#}
#################################################################################
#
#
#
#### Get protein - lipid overlaps ###############################################
#my @resOverlapScore;
#my @resCavityScore;
#my @resOverlapList;
#my @resCavityList;
#my @gridSlices;
#
#foreach (@inflGroupIds) {
#    my $zSlice = 0.5; # A z-slice with a width of 0.5 nm give good results.
#
#    for (my $zSliceMin=$lipidData[$_]{'zMin'}; $zSliceMin<=$lipidData[$_]{'zMax'}; $zSliceMin+=$zSlice) {
#        ### Get the protein grid of the current z-slice ########################
#        my $zSliceMax = $zSliceMin + $zSlice - 0.001; # Three numbers behind the point (GRO file accuracy).
#        my $zSliceStr = $zSliceMin . "-" . $zSliceMax;
#        print "    Running through the z-slice (range $zSliceStr)\n";
#
#        my @tmpArray = Protein::analyzeProtein(\@ndxData,
#                                               $protGroupId,
#                                               \@suGroupIds,
#                                               $coordData{'atoms'},
#                                               $coordData{'box'},
#                                               $gridDelta,
#                                               $zSliceMin,
#                                               $zSliceMax,
#                                               $extCavDetection,
#                                               $ndxData[$_]{'groupName'} . ".zslice_$zSliceStr.xy");
#        $gridSlices[$_]{$zSliceStr} = $tmpArray[2];
#        ########################################################################
#
#        getOverlaps($gridSlices[$_]{$zSliceStr},
#                    $ndxData[$_]{'atoms'},
#                    $coordData{'atoms'},
#                    $zSliceMax,
#                    $zSliceMin,
#                    \@{$resOverlapScore[$_]},
#                    \@{$resCavityScore[$_]},
#                    $lipidFilledCav);
#    }
##    overlapScore2PdbBeta("ols2beta.pdb", \%coordData, $resOverlapScore[$_]);
#
#    my @resOverlapScoreSorted = sortResOverlaps($resOverlapScore[$_]);
#    my @resCavityScoreSorted  = sortResOverlaps($resCavityScore[$_]);
#
#
#    ### Part the detected lipids to their leaflets #############################
#    foreach my $item (@resOverlapScoreSorted) {
#        if ($lipidData[$_]{'leaflets'}{'upper'}{'uResIds'}[ $$item{'uResId'} ]) {
#            push(@{$resOverlapList[$_]{'upper'}}, $$item{'uResId'});
#        }
#        elsif ($lipidData[$_]{'leaflets'}{'lower'}{'uResIds'}[ $$item{'uResId'} ]) {
#            push(@{$resOverlapList[$_]{'lower'}}, $$item{'uResId'});
#        }
#    }
#
#    foreach my $item (@resCavityScoreSorted) { # Cavity lipids.
#        if ($lipidData[$_]{'leaflets'}{'upper'}{'uResIds'}[ $$item{'uResId'} ]) {
#            push(@{$resCavityList[$_]{'upper'}}, $$item{'uResId'});
#        }
#        elsif ($lipidData[$_]{'leaflets'}{'lower'}{'uResIds'}[ $$item{'uResId'} ]) {
#            push(@{$resCavityList[$_]{'lower'}}, $$item{'uResId'});
#        }
#    }
#    ############################################################################
#}
#################################################################################
#
#
#
#my @uResIdAtomIdList = uResId2AtomIdList($coordData{'atoms'}); # Get the atom IDs of the overlapping residues.
#my %protGeoCenter;
#my @resGeoCenter;
#my @inflAtomList; # Dynamic components of each group (lipids for inflating).
#my @statAtomList; # Static components of each group.
#my @cavAtomList;  # Cavity components of each group.
#
#my @nLipidCav;
#my @nLipidDel;
#
#my %statNdxGroup = ('groupName' => "Static");
#my %dynaNdxGroup = ('groupName' => "Dynamic");
#my %cavNdxGroup  = ('groupName' => "Cavity");
#
#%protGeoCenter = getGeoCenterXY($ndxData[$protGroupId]{'atoms'}, $coordData{'atoms'});
#
#foreach (@inflGroupIds) {
#    my %cavRes;
#    my %delRes;
#
#    ### Cavity lipid part ######################################################
#    if ($lipidFilledCav) {
#        $nLipidCav[$_]{'upper'} = getnLipPerArea($gridACav[$_]{'upper'}, $coordData{'box'}, $lipidData[$_]{'leaflets'}{'upper'}{'nLipids'});
#        $nLipidCav[$_]{'lower'} = getnLipPerArea($gridACav[$_]{'lower'}, $coordData{'box'}, $lipidData[$_]{'leaflets'}{'lower'}{'nLipids'});
#        print "Upper leaflet: $nLipidCav[$_]{'upper'} ($gridAProt[$_]{'upper'}); Lower leaflet: $nLipidCav[$_]{'lower'} ($gridAProt[$_]{'lower'})\n";
#        print "Upper leaflet: " . int($nLipidCav[$_]{'upper'}) . "; Lower leaflet: " . int($nLipidCav[$_]{'lower'}) . "\n";
#
#        print "Cavity lipids in the upper leaflet:\n";
#        for (my $i=0; $i<int($nLipidCav[$_]{'upper'}); $i++) {
#            print " " . $resCavityList[$_]{'upper'}[$i];
#            $cavRes{'upper'}[$resCavityList[$_]{'upper'}[$i]] = 1;
#        }
#        print "\n";
#
#        print "Cavity lipids in the lower leaflet:\n";
#        for (my $i=0; $i<int($nLipidCav[$_]{'lower'}); $i++) {
#            print " " . $resCavityList[$_]{'lower'}[$i];
#            $cavRes{'lower'}[$resCavityList[$_]{'lower'}[$i]] = 1;
#        }
#        print "\n";
#    }
#    ############################################################################
#
#
#    ### Delete residues ########################################################
#    $nLipidDel[$_]{'upper'} = getnLipPerArea($gridAProt[$_]{'upper'}, $coordData{'box'}, $lipidData[$_]{'leaflets'}{'upper'}{'nLipids'});
#    $nLipidDel[$_]{'lower'} = getnLipPerArea($gridAProt[$_]{'lower'}, $coordData{'box'}, $lipidData[$_]{'leaflets'}{'lower'}{'nLipids'});
#    print "Upper leaflet: $nLipidDel[$_]{'upper'} ($gridAProt[$_]{'upper'}); Lower leaflet: $nLipidDel[$_]{'lower'} ($gridAProt[$_]{'lower'})\n";
#    print "Upper leaflet: " . round($nLipidDel[$_]{'upper'}, 1) . "; Lower leaflet: " . round($nLipidDel[$_]{'lower'}, 1) . "\n";
#
#    print "Will delete lipids in the upper leaflet:\n";
#    my $delOffset = 0;
#    for (my $i=0; $i<round($nLipidDel[$_]{'upper'}, 1); $i++) {
#        my $iTmp = $i+$delOffset;
#        while ($cavRes{'upper'}[$resOverlapList[$_]{'upper'}[$iTmp]]) {
#            $delOffset++;
#            $iTmp = $i+$delOffset;
#        }
#        print " " . $resOverlapList[$_]{'upper'}[$iTmp];
#        $delRes{'upper'}[$resOverlapList[$_]{'upper'}[$iTmp]] = 1;
#    }
#    print "\n";
#
#    print "Will delete lipids in the lower leaflet:\n";
#    $delOffset = 0;
#    for (my $i=0; $i<round($nLipidDel[$_]{'lower'}, 1); $i++) {
#        my $iTmp = $i+$delOffset;
#        while ($cavRes{'lower'}[$resOverlapList[$_]{'lower'}[$iTmp]]) {
#            $delOffset++;
#            $iTmp = $i+$delOffset;
#        }
#        print " " . $resOverlapList[$_]{'lower'}[$iTmp];
#        $delRes{'lower'}[$resOverlapList[$_]{'lower'}[$iTmp]] = 1;
#    }
#    print "\n";
#
#
#    my @tmpArray = delResidues($coordData{'atoms'}, $delRes{'upper'});
#    $coordData{'atoms'} = $tmpArray[0];
#    my $delAtomsDataRef = $tmpArray[1];
#    updateDelResnames(\%delResnames, $delAtomsDataRef) if $topFile && !$noDeflating;
#
#    @tmpArray = delResidues($coordData{'atoms'}, $delRes{'lower'});
#    $coordData{'atoms'} = $tmpArray[0];
#    $delAtomsDataRef    = $tmpArray[1];
#    updateDelResnames(\%delResnames, $delAtomsDataRef) if $topFile && !$noDeflating;
#    ############################################################################
#
#
#    my @inflResList; # Only to check the output. COMMENT THIS OUT.
#    my @statResList; # Only to check the output. COMMENT THIS OUT.
#    my @cavResList;  # Only to check the output. COMMENT THIS OUT.
#    for (my $uResId=0; $uResId<=@{$resOverlapScore[$_]}; $uResId++) {
#        next unless $uResIdAtomIdList[$uResId][0];
#        next unless $coordData{'atoms'}[$uResIdAtomIdList[$uResId][0]];
#        next unless defined $resOverlapScore[$_][$uResId];
#        unless ($resOverlapScore[$_][$uResId]) {
#            push(@{$statAtomList[$_]}, @{$uResIdAtomIdList[$uResId]});
#            push(@statResList, $uResId); # COMMENT THIS OUT.
#            next;
#        }
#
#        ### Get the geometrical center of each inflatable lipid ################
#        %{$resGeoCenter[$uResId]} = getGeoCenterXY($uResIdAtomIdList[$uResId], $coordData{'atoms'});
##        printf("Geometrical center of the reference lipid %d is: x = %.2f, y=%.2f\n", $resId, $resGeoCenter[$resId]{'cooX'}, $resGeoCenter[$resId]{'cooY'});
#        ########################################################################
#
#        if ($lipidFilledCav && ($cavRes{'upper'}[$uResId] || $cavRes{'lower'}[$uResId])) {
#            push(@{$cavAtomList[$_]}, @{$uResIdAtomIdList[$uResId]});
#            push(@cavResList, $uResId); # Only to check the output. COMMENT THIS OUT.
#
#            if ($cavRes{'upper'}[$uResId]) {
#                foreach my $cavAtomId (@{$uResIdAtomIdList[$uResId]}) {
#                    $groupCavAtomInflVectors[$_][$cavAtomId]{'cooX'} = $coordData{'atoms'}[$cavAtomId]{'cooX'} - $cavGeoCenterXy[$_]{'upper'}{'cooX'};
#                    $groupCavAtomInflVectors[$_][$cavAtomId]{'cooY'} = $coordData{'atoms'}[$cavAtomId]{'cooY'} - $cavGeoCenterXy[$_]{'upper'}{'cooY'};
#                    $groupCavAtomInflVectors[$_][$cavAtomId]{'cooZ'} = $coordData{'atoms'}[$cavAtomId]{'cooZ'} - $cavGeoCenterXy[$_]{'upper'}{'cooZ'};
#
#                    $coordData{'atoms'}[$cavAtomId]{'cooX'} = $cavGeoCenterXy[$_]{'upper'}{'cooX'};
#                    $coordData{'atoms'}[$cavAtomId]{'cooY'} = $cavGeoCenterXy[$_]{'upper'}{'cooY'};
#                    $coordData{'atoms'}[$cavAtomId]{'cooZ'} = $cavGeoCenterXy[$_]{'upper'}{'cooZ'};
#                }
#            }
#            elsif ($cavRes{'lower'}[$uResId]) {
#                foreach my $cavAtomId (@{$uResIdAtomIdList[$uResId]}) {
#                    print "Ascribe atom $cavAtomId a translation vector\n";
#                    $groupCavAtomInflVectors[$_][$cavAtomId]{'cooX'} = $coordData{'atoms'}[$cavAtomId]{'cooX'} - $cavGeoCenterXy[$_]{'lower'}{'cooX'};
#                    $groupCavAtomInflVectors[$_][$cavAtomId]{'cooY'} = $coordData{'atoms'}[$cavAtomId]{'cooY'} - $cavGeoCenterXy[$_]{'lower'}{'cooY'};
#                    $groupCavAtomInflVectors[$_][$cavAtomId]{'cooZ'} = $coordData{'atoms'}[$cavAtomId]{'cooZ'} - $cavGeoCenterXy[$_]{'lower'}{'cooZ'};
#
#                    $coordData{'atoms'}[$cavAtomId]{'cooX'} = $cavGeoCenterXy[$_]{'lower'}{'cooX'};
#                    $coordData{'atoms'}[$cavAtomId]{'cooY'} = $cavGeoCenterXy[$_]{'lower'}{'cooY'};
#                    $coordData{'atoms'}[$cavAtomId]{'cooZ'} = $cavGeoCenterXy[$_]{'lower'}{'cooZ'};
#                }
#            }
#        }
#        else {
#            push(@{$inflAtomList[$_]}, @{$uResIdAtomIdList[$uResId]});
#            push(@inflResList, $uResId); # Only to check the output. COMMENT THIS OUT.
#        }
#    }
#    push(@{$statNdxGroup{'atoms'}}, @{$statAtomList[$_]});
#    push(@{$dynaNdxGroup{'atoms'}}, @{$inflAtomList[$_]});
#    push(@{$cavNdxGroup{'atoms'}}, @{$cavAtomList[$_]}) if $lipidFilledCav;
#
#    print "Will inflate these residues\n";
#    for (my $i=0; $i<@inflResList; $i++) {
#        print " " . $inflResList[$i];
#    }
#    print "\n";
#    print "Will handle these residues as static\n";
#    for (my $i=0; $i<@statResList; $i++) {
#        print " " . $statResList[$i];
#    }
#    print "\n";
#
#    if ($lipidFilledCav) {
#        print "Will handle these residues as cavity\n";
#        for (my $i=0; $i<@cavResList; $i++) {
#            print " " . $cavResList[$i];
#        }
#        print "\n";
#    }
#    
#    ### Recursive function for dynamic inflating ###############################
#    $groupInflatingFactors[$_] = dynInflating($gridSlices[$_], $inflAtomList[$_], $coordData{'atoms'}, \%protGeoCenter, \@resGeoCenter);
#    print "Used an inflating factor of " . $groupInflatingFactors[$_] . " to solve protein - lipid overlaps\n";
#    ############################################################################
#}
#
#
#push(@ndxData, \%statNdxGroup);
#push(@ndxData, \%dynaNdxGroup);
#push(@ndxData, \%cavNdxGroup) if $lipidFilledCav;
#
#
#
#### Update internal data after lipid deletion ##################################
#my @newNdxData;
#my @renumMatrix;
#($coordData{'atoms'}, @renumMatrix) = renumAtoms($coordData{'atoms'});
#$coordData{'nAtoms'} = scalar(@{$coordData{'atoms'}}) - 1;
#
##unless ($lipidResname && $noDeflating) {
#    @newNdxData = NDXFiles::updateNdxData(\@ndxData, \@renumMatrix);
##    NDXFiles::printNdxGroups(@newNdxData); # Just check the updated NDX data.
#    NDXFiles::writeNdx("inflated.ndx", \@newNdxData);
##}
#################################################################################
#
#
#
#### Inflate the box ############################################################
#my $maxInflFactor = 0;
#foreach (@inflGroupIds) {
#    $maxInflFactor = $groupInflatingFactors[$_] if $groupInflatingFactors[$_] > $maxInflFactor;
#}
#inflateBoxXY($coordData{'box'}, $maxInflFactor);
#centerGroupXY($newNdxData[0]{'atoms'}, $coordData{'atoms'}, $coordData{'box'});
#################################################################################
#
#
#
#### Write inflated coordinates to an output coordinates file ###################
#if ($outFile =~ /\.pdb$/) {
#    PDBFiles::writePdb($outFile, \%coordData);
#}
#else {
#    GROFiles::writeGro($outFile, \%coordData);
#}
#################################################################################
#
#
#
#### Update topology ############################################################
#if ($topFile) {
#    %topData = TOPFiles::readTop($topFile); # Read input TOP file.
#    foreach my $resname (keys %delResnames) {
#        my $tmpDir = '.';
#        $tmpDir = $1 if ($topFile =~ /^(.*)\/[^\/]+$/);
#        my $molType = resname2moltype($topData{'include'}, $resname, $tmpDir . "/");
#        TOPFiles::updateTop($topFile, $molType, $topData{'molecules'}{$resname} - $delResnames{$resname}) if $molType;
#    }
#}
#################################################################################
#
#
#
#### Shrink the model ###########################################################
#exit if $noDeflating;
#
#
#### Get the deflating factor of each group #####################################
#foreach (@inflGroupIds) {
#    $groupDeflatingFactors[$_] = getDeflatingFactor($groupInflatingFactors[$_], $nDeflatingSteps);
#    print "Will use a deflating factor of " . $groupDeflatingFactors[$_] . " to shrink the lipids back in $nDeflatingSteps steps\n";
#}
#################################################################################
#
#
#print "Deflating\n";
#
#my @newInflAtomList;
#my @newCavAtomList;
#foreach my $groupId (@inflGroupIds) {
#    foreach (@{$inflAtomList[$groupId]}) {
#        push(@{$newInflAtomList[$groupId]}, $renumMatrix[$_]);
#    }
#    foreach (@{$cavAtomList[$groupId]}) {
#        push(@{$newCavAtomList[$groupId]}, $renumMatrix[$_]);
#        print "Will convert atom Id $_ to $renumMatrix[$_]\n";
#    }
#}
#
#
#my $shrinkOutFile = sprintf("shrink.%02d.gro", 0);
#doEM($mdpFile, $outFile, $topFile, "inflated.ndx", $shrinkOutFile);
##GROFiles::writeGro($shrinkOutFile, \%coordData);
#
#for (my $step=0; $step<$nDeflatingSteps; $step++) {
#    $shrinkOutFile = sprintf("shrink.%02d.gro", $step);
#    %coordData = GROFiles::readGro($shrinkOutFile); # Read input GRO file.
#
#    ### Dynamic box sice handling (multiple bilayers) ##########################
#    my %tmpBoxXY = ('cooX' => 0, 'cooY' => 0);
#    foreach (@inflGroupIds) {
#        my $boxX = $coordData{'box'}{'cooX'} * $groupDeflatingFactors[$_];
#        my $boxY = $coordData{'box'}{'cooY'} * $groupDeflatingFactors[$_];
#        $tmpBoxXY{'cooX'} = $boxX if $boxX > $tmpBoxXY{'cooX'};
#        $tmpBoxXY{'cooY'} = $boxY if $boxY > $tmpBoxXY{'cooY'};
#    }
#    $coordData{'box'}{'cooX'} = $tmpBoxXY{'cooX'};
#    $coordData{'box'}{'cooY'} = $tmpBoxXY{'cooY'};
#    ############################################################################
#
#    %protGeoCenter = getGeoCenterXY($newNdxData[$protGroupId]{'atoms'}, $coordData{'atoms'});
#
#    foreach (@inflGroupIds) {
#        my @resGeoCenter = getInflGeoCenterXY($newInflAtomList[$_], $coordData{'atoms'});
#        scaleGroupXY($newInflAtomList[$_], $coordData{'atoms'}, $groupDeflatingFactors[$_], \%protGeoCenter, \@resGeoCenter, 1);
#
#        if ($lipidFilledCav) {
#            my $tmpFactor = 1/$nDeflatingSteps;
#            foreach my $cavAtomId (@{$cavAtomList[$_]}) {
#                print "Will inflate cavity lipid atom $cavAtomId\n";
#                $coordData{'atoms'}[ $renumMatrix[$cavAtomId] ]{'cooX'} += $groupCavAtomInflVectors[$_][$cavAtomId]{'cooX'} * $tmpFactor;
#                $coordData{'atoms'}[ $renumMatrix[$cavAtomId] ]{'cooY'} += $groupCavAtomInflVectors[$_][$cavAtomId]{'cooY'} * $tmpFactor;
#                $coordData{'atoms'}[ $renumMatrix[$cavAtomId] ]{'cooZ'} += $groupCavAtomInflVectors[$_][$cavAtomId]{'cooZ'} * $tmpFactor;
#            }
#        }
#    }
#    centerGroupXY($newNdxData[0]{'atoms'}, $coordData{'atoms'}, $coordData{'box'});
#
#    $shrinkOutFile = sprintf("shrink.%02d.gro", $step+1);
##    GROFiles::writeGro($shrinkOutFile, \%coordData);
#    GROFiles::writeGro('tmp_out.gro', \%coordData);
#    doEM($mdpFile, 'tmp_out.gro', $topFile, "inflated.ndx", $shrinkOutFile);
#}


### Update NDX file for the last energy minimization step ######################
#pop(@newNdxData);
#pop(@newNdxData);
#my @empty;
#my @dynaAtomList;
#foreach (@statAtomList) {
#    push(@dynaAtomList, $renumMatrix[$_]);
#}
#foreach (@inflAtomList) {
#    push(@dynaAtomList, $renumMatrix[$_]);
#}
#my %statNdxGroup = ('groupName' => "Static");
#my %dynaNdxGroup = ('groupName' => "Dynamic");
#$statNdxGroup{'atoms'}[0] = 1; # NOTE: Atom number one is set as static! This may lead to problems, but the group must not be empty.
#$dynaNdxGroup{'atoms'} = \@dynaAtomList;
#push(@newNdxData, \%statNdxGroup);
#push(@newNdxData, \%dynaNdxGroup);
#my @embeddedNdxData = NDXFiles::updateNdxData(\@newNdxData, \@renumMatrix);
#NDXFiles::printNdxGroups(@embeddedNdxData); # Just check the updated NDX data.
#NDXFiles::writeNdx("embedded.ndx", \@embeddedNdxData);
################################################################################

#doEM($mdpFile, $shrinkOutFile, $topFile, "embedded.ndx", "embedded.gro");
################################################################################


sub getInflGeoCenterXY {
    my $inflAtomListRef = shift;
    my $coordsRef       = shift;
    my @resGeoCenter;
    my @residAtomCounter;
    my @resIdCounter;
    my @uResIdList;

    foreach (@{$inflAtomListRef}) {
        next unless $$coordsRef[$_];
        my $uResId = $$coordsRef[$_]{'uResId'};
        $resGeoCenter[$uResId]{'cooX'} += $$coordsRef[$_]{'cooX'};
        $resGeoCenter[$uResId]{'cooY'} += $$coordsRef[$_]{'cooY'};
        $residAtomCounter[$uResId]++;
        next if $resIdCounter[$uResId];
        push(@uResIdList, $uResId);
        $resIdCounter[$uResId] = 1;
    }
    foreach (@uResIdList) {
        $resGeoCenter[$_]{'cooX'} /= $residAtomCounter[$_];
        $resGeoCenter[$_]{'cooY'} /= $residAtomCounter[$_];
    }
    return @resGeoCenter;
}


printFoot();
exit;


sub overlapScore2PdbBeta {
    my $pdbFile            = shift;
    my $coordDataRef       = shift;
    my $resOverlapScoreRef = shift;

    foreach (@{$$coordDataRef{'atoms'}}) {
#        print $$resOverlapScoreRef[ $$_{'uResId'} ] . " $$_{'uResId'}\n";
        $$_{'tempFactor'} = $$resOverlapScoreRef[ $$_{'uResId'} ];
    }

    PDBFiles::writePdb($pdbFile, $coordDataRef);
}



sub getDeflatingFactor {
    my $inflatingFactor = shift;
    my $nDeflatingSteps = shift;

    return exp((log(1/$inflatingFactor))/$nDeflatingSteps);
}



sub dynInflating {
    my $gridSlicesRef    = shift;
    my $inflAtomListRef  = shift;
    my $coordsRef        = shift;
    my $protGeocenterRef = shift;
    my $resGeoCenterRef  = shift;
    my $scalingFactor    = shift;

    $scalingFactor = 1.1 unless $scalingFactor;


    ### Inflate the defined lipids #############################################
    my $scaledCoordsRef = scaleGroupXY($inflAtomListRef, $coordsRef, $scalingFactor, $protGeocenterRef, $resGeoCenterRef);
    ############################################################################


    ### Check the inflated lipid atoms for protein overlaps ####################
    unless (checkOverlap($gridSlicesRef, $inflAtomListRef, $scaledCoordsRef)) {
        foreach (@{$inflAtomListRef}) {
            $$coordsRef[$_]{'cooX'} = $$scaledCoordsRef[$_]{'cooX'};
            $$coordsRef[$_]{'cooY'} = $$scaledCoordsRef[$_]{'cooY'};
        }
        return sprintf("%.1f", $scalingFactor);
    }
    ############################################################################


    ### Still overlaps? Inflate again! #########################################
    return dynInflating($gridSlicesRef, $inflAtomListRef, $coordsRef, $protGeocenterRef, $resGeoCenterRef, $scalingFactor+0.1);
    ############################################################################
}



sub scaleGroupXY {
    my $atomIdsRef       = shift;
    my $coordsRef        = shift;
    my $scalingFactor    = shift;
    my $protGeoCenterRef = shift;
    my $resGeoCenterRef  = shift;
    my $direct           = shift; # Apply translation directly to the coordinates (1) or create new coordinates (0)?
    my @translVec; # Scaling vector (x,y) for each lipid residue.
    my @scaledCoords;


    foreach (@{$atomIdsRef}) {
        my $uResId = $$coordsRef[$_]{'uResId'};
        next unless $$resGeoCenterRef[$uResId];
        next if $translVec[$uResId];
        $translVec[$uResId]{'cooX'} = ($$resGeoCenterRef[$uResId]{'cooX'} - $$protGeoCenterRef{'cooX'}) * $scalingFactor - ($$resGeoCenterRef[$uResId]{'cooX'} - $$protGeoCenterRef{'cooX'});
        $translVec[$uResId]{'cooY'} = ($$resGeoCenterRef[$uResId]{'cooY'} - $$protGeoCenterRef{'cooY'}) * $scalingFactor - ($$resGeoCenterRef[$uResId]{'cooY'} - $$protGeoCenterRef{'cooY'});
    }


    foreach (@{$atomIdsRef}) {
        my $uResId = $$coordsRef[$_]{'uResId'};
        next unless $translVec[$uResId]{'cooX'};
        if ($direct) {
            $$coordsRef[$_]{'cooX'} += $translVec[$uResId]{'cooX'};
            $$coordsRef[$_]{'cooY'} += $translVec[$uResId]{'cooY'};
        }
        else {
            $scaledCoords[$_]{'atomName'} = $$coordsRef[$_]{'atomName'};
            $scaledCoords[$_]{'cooX'} = $$coordsRef[$_]{'cooX'} + $translVec[$uResId]{'cooX'};
            $scaledCoords[$_]{'cooY'} = $$coordsRef[$_]{'cooY'} + $translVec[$uResId]{'cooY'};
            $scaledCoords[$_]{'cooZ'} = $$coordsRef[$_]{'cooZ'};
        }
    }
    return \@scaledCoords;
}



#sub checkOverlap {
#    my $gridSlicesRef = shift;
#    my $atomsIdsRef   = shift;
#    my $coordsRef     = shift;
#
#    foreach my $zSliceStr (keys %{$gridSlicesRef}) {
#        next unless $zSliceStr =~ /^([-+]?\d*\.?\d+([eE][-+]?\d+)?)-([-+]?\d*\.?\d+([eE][-+]?\d+)?)$/;
#        my $zSliceMin = $1;
#        my $zSliceMax = $3;
#
#        my $gridsizeX = @{$$gridSlicesRef{$zSliceStr}};
#        my $gridsizeY = @{$$gridSlicesRef{$zSliceStr}[0]};
#
#        foreach (@{$atomsIdsRef}) {
#            next unless $$coordsRef[$_]{'cooZ'};
#            next if $$coordsRef[$_]{'cooZ'} < $zSliceMin;
#            next if $$coordsRef[$_]{'cooZ'} > $zSliceMax;
#
#            my $element  = substr($$coordsRef[$_]{'atomName'}, 0, 1);
#            my $radius   = PTE::getRadius($element);
#            my $radiusSq = $radius * $radius;
#
#            my $tmpGridX = int($$coordsRef[$_]{'cooX'} / $gridDelta);
#            my $tmpGridY = int($$coordsRef[$_]{'cooY'} / $gridDelta);
#            my $subrange = int($radiusSq / $gridDelta);
#
#            for (my $x=getMax($tmpGridX-$subrange, 0); $x<=getMin($tmpGridX+$subrange, $gridsizeX-1); $x++) {
#                for (my $y=getMax($tmpGridY-$subrange, 0); $y<=getMin($tmpGridY+$subrange, $gridsizeY-1); $y++) {
#                    next unless $$gridSlicesRef{$zSliceStr}[$x][$y];
#
#                    my $dx = $$coordsRef[$_]{'cooX'} - $x * $gridDelta;
#                    my $dy = $$coordsRef[$_]{'cooY'} - $y * $gridDelta;
#                    next if ($dx*$dx + $dy*$dy) > $radiusSq;
#
#                    return 1 if ($$gridSlicesRef{$zSliceStr}[$x][$y]);
#                }
#            }
#        }
#    }
#    return 0;
#}
#
#
#
#sub delResidues {
#    my $coordsRef = shift;
#    my $delResRef = shift;
#    my @restAtomsArray;
#    my @delAtomsArray;
#
#    for (my $i=0; $i<@{$coordsRef}; $i++) {
#        next unless $$coordsRef[$i]{'uResId'};
##        print "Will remove residue $$coordsRef[$i]{'uResId'} (atom $i)\n";
#        $$delResRef[$$coordsRef[$i]{'uResId'}] ? $delAtomsArray[$i] = $$coordsRef[$i] : $restAtomsArray[$i] = $$coordsRef[$i];
#    }
#    return (\@restAtomsArray, \@delAtomsArray);
#}
#
#
#
#sub uResId2AtomIdList {
#    my $coordsRef = shift;
#
#    my @uResIdAtomIdList;
#
#    for (my $atomId=0; $atomId<@{$coordsRef}; $atomId++) {
#        next unless $$coordsRef[$atomId]{'uResId'};
#        push(@{$uResIdAtomIdList[$$coordsRef[$atomId]{'uResId'}]}, $atomId);
#    }
#    return @uResIdAtomIdList;
#}
#
#
#
#sub getnLipPerArea {
#    my $areaProtein = shift;
#    my $boxRef      = shift;
#    my $nLipids     = shift;
#
#    my $areaPerLipid = $$boxRef{'cooX'} * $$boxRef{'cooY'} / $nLipids;
#
#    return ($areaProtein / $areaPerLipid);
#}
#
#
#
#sub sortResOverlaps {
#    my $resOverlapScoreRef = shift;
#
#    my @tmp;
#    for (my $resId=0; $resId<@{$resOverlapScoreRef}; $resId++) {
#        next unless $$resOverlapScoreRef[$resId];
#        my %tmpHash = ('uResId' => $resId,
#                       'score'  => $$resOverlapScoreRef[$resId]);
#        push(@tmp, \%tmpHash);
#    }
#    return sort({$b->{'score'}<=>$a->{'score'}} @tmp);
#}
#
#
#
#sub getOverlaps {
#    my $gridRef            = shift;
#    my $atomsIdsRef        = shift;
#    my $coordsRef          = shift;
#    my $zLimitUpper        = shift;
#    my $zLimitLower        = shift;
#    my $resOverlapScoreRef = shift;
#    my $resCavityScoreRef  = shift;
#    my $lipidFilledCav     = shift;
#
#    my $gridsizeX = @$gridRef;
#    my $gridsizeY = @$gridRef[0];
#    my $nOverlaps   = 0;
#    my $atomCounter = 0;
#    my @atomsWithOverlapScore;
#
#    my @specialGrid;
#
#    print "      Searching for protein - lipid overlaps: 0.000%\r" if $main::verbose;
#    foreach (@{$atomsIdsRef}) {
#        next unless $$coordsRef[$_]{'cooZ'};
#        printf("      Searching for protein - lipid overlaps: %d%% (found overlaps %d)\r", ($atomCounter++)/scalar(@{$atomsIdsRef})*100, $nOverlaps) if $main::verbose;
#        next if $$coordsRef[$_]{'cooZ'} > $zLimitUpper;
#        next if $$coordsRef[$_]{'cooZ'} < $zLimitLower;
#
#        my $uResId = $$coordsRef[$_]{'uResId'};
#        $$resOverlapScoreRef[$uResId] = 0 unless $$resOverlapScoreRef[$uResId];
#        $$resCavityScoreRef[$uResId]  = 0 unless $$resCavityScoreRef[$uResId];
#
#        my $element = substr($$coordsRef[$_]{'atomName'}, 0, 1);
#        my $radius  = PTE::getRadius($element);
#        my $radius2 = $radius * $radius;
#
#        my $tmpGridX = int($$coordsRef[$_]{'cooX'} / $gridDelta);
#        my $tmpGridY = int($$coordsRef[$_]{'cooY'} / $gridDelta);
#        my $subrange = int($radius / $gridDelta);
#
#        for (my $x=getMax($tmpGridX-$subrange, 0); $x<=getMin($tmpGridX+$subrange, $gridsizeX-1); $x++) {
#            for (my $y=getMax($tmpGridY-$subrange, 0); $y<=getMin($tmpGridY+$subrange, $gridsizeY-1); $y++) {
#                next unless $$gridRef[$x][$y];
#
#                my $dx = $$coordsRef[$_]{'cooX'} - $x * $gridDelta;
#                my $dy = $$coordsRef[$_]{'cooY'} - $y * $gridDelta;
#                next if ($dx*$dx + $dy*$dy) > $radius2;
#
#                if ($$gridRef[$x][$y] == 2) {
#                    $$resOverlapScoreRef[$uResId] += $lipidFilledCav ? 1 : 4; # Prefer cavity lipids for deletion if there is no lipid-filled cavity.
#                    $$resCavityScoreRef[$uResId]++;
#                }
#                else {
#                    $$resOverlapScoreRef[$uResId]++;
#                }
#
#                $nOverlaps++ unless $atomsWithOverlapScore[$_];
#                $atomsWithOverlapScore[$_] = 1;
#            }
#        }
#    }
#    printf("      Searching for protein - lipid overlaps: 100%% (found overlaps %d)\n", $nOverlaps) if $main::verbose;
#}



#sub doEM {
#    my $mdpFile   = shift;
#    my $coordFile = shift;
#    my $topFile   = shift;
#    my $ndxFile   = shift;
#    my $shrinkOut = shift;
#
#    my $tmpCmd = sprintf("grompp -f %s -c %s -p %s -n %s -o %s -maxwarn 1", $mdpFile, $coordFile, $topFile, $ndxFile, 'tmp.tpr');
#    print $tmpCmd . "\n";
#    `$tmpCmd 1>>log.out 2>&1`;
#
#    $tmpCmd = sprintf("mdrun -s %s -v -deffnm %s -c %s", 'tmp.tpr', 'tmp_out', $shrinkOut);
#    print $tmpCmd . "\n";
#    `$tmpCmd 1>>log.out 2>&1`;
#
#    system("rm tmp.tpr tmp_out.* step*.pdb mdout.mdp");
#}



################################################################################
### Subroutines ################################################################
################################################################################
sub printHead {
    my @headLines = ("################################################################################",
                     "",
                     "Tango $version",
                     "Molecular thickness measurement using atomic van der Waals radii.",
                     "Copyright Thomas H. Schmidt, $year",
                     "",
                     "http://code.google.com/p/dancesteps",
                     "",
                     "Tango comes with ABSOLUTELY NO WARRANTY.",
                     "This is free software, and you are welcome to redistribute it",
                     "under certain conditions; type `-copyright' for details.",
                     "",
                     "################################################################################");
    my $maxLength = 80;
    foreach (@headLines) {
        $maxLength = (length $_ > $maxLength) ? length($_) : $maxLength;
    }

    foreach (@headLines) {
        printf "%s%-${maxLength}s\n", ' ' x int(($maxLength - length($_))/2), $_;
    }
}



sub printFoot {
    print <<EndOfFoot;
Please cite:
  [1] Schmidt, T. H. Tango: Grid-based Membrane Analysis (Manual)

EndOfFoot
}



sub printHelp {
    my $cmdLParamRef   = shift;
    my $quitAfterPrint = shift;


    print <<EndOfHelp;
DESCRIPTION
-----------
Tango reads a GRO or PDB coordinate file and calculates the diameter of a group
defined by a given GROMACS NDX file or by residue name.
Therefor Tango maps the system to a 3D grid of a given resolution (-dx, -dy, -dz).

USAGE: tango -f GROFILE -o OUTPUTFILE

EndOfHelp

    printParamHelp($cmdLParamRef);

    printFoot();

    exit if $quitAfterPrint;
}



sub printCopyright {
    print <<"EndOfCopyright";
Tango is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
any later version.

Tango is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with Tango; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

EndOfCopyright
    exit;
}



#sub cmdlineParser {
#    my $paramsRef = shift;
#
#    my %knownParam;
#    my @unknownParam;
#
#    for (my $argID=0; $argID<@ARGV; $argID++) {
#        my $cmdlineName = $ARGV[$argID];
#        $knownParam{$ARGV[$argID]} = 0;
#
#        foreach my $paramKey (keys %{$paramsRef}) {
#            my $paramName = $paramKey;
#            my $paramType = 0;
#            my $paramIni  = "--";
#            my $isArray   = 0;
#
#            if ($paramKey =~ /^(.+?)=([\@f])$/) {
#                $paramName = $1;
#                $paramType = $2;
#            }
#
#            $paramIni = "-" if (length($paramName) == 1);
#
#            if ($paramType eq "@") {
#                $isArray = 1;
#            }
#            elsif ($paramType eq "f") {
#                if ($cmdlineName eq $paramIni.$paramName) {
#                    if (ref(${$paramsRef}{$paramKey}) eq "SCALAR") {
#                        ${${$paramsRef}{$paramKey}} = 1;
#                        $knownParam{$cmdlineName} = 1;
#                    }
#                    elsif (ref(${$paramsRef}{$paramKey}) eq "CODE") {
#                        &{${$paramsRef}{$paramKey}}();
#                        $knownParam{$cmdlineName} = 1;
#                    }
#                }
#                next;
#            }
#
#            if ($cmdlineName eq $paramIni.$paramName && not $isArray) {
#                $argID++;
#                next if($ARGV[$argID] =~ /^-/);
#
#                ${${$paramsRef}{$paramKey}} = $ARGV[$argID];
#                $knownParam{$cmdlineName} = 1;
#            }
#            elsif ($ARGV[$argID] eq $paramIni.$paramName && $isArray) {
#                $knownParam{$cmdlineName} = 1;
#                $argID++;
#
#                my @tmpArray;
#                while ($argID <= $#ARGV && not $ARGV[$argID] =~ /^--/ && not $ARGV[$argID] =~ /^-.$/) {
#                    push (@tmpArray, $ARGV[$argID]);
#                    $argID++;
#                }
#                @{${$paramsRef}{$paramKey}} = @tmpArray;
#                $argID--;
#            }
#        }
#
#        if (defined $knownParam{$cmdlineName} && $knownParam{$cmdlineName} == 0) {
#            push(@unknownParam, $cmdlineName) if($cmdlineName =~ /^-/);
#        }
#    }
#
#
#    ### Catch unknown parameters ###############################################
#    if (@unknownParam && ${$paramsRef}{"UNKNOWN"} && ref(${$paramsRef}{"UNKNOWN"}) eq "CODE") {
#        print "WARNING: Unknown or non-set parameters detected:\n";
#        for (@unknownParam) { print "        \"$_\"\n"; }
#        &{${$paramsRef}{"UNKNOWN"}}();
#    }
#    elsif (@unknownParam) {
#        print "ERROR: Unknown or non-set parameters detected:\n";
#        for(@unknownParam) { print "       \"$_\"\n"; }
#        exit;
#    }
#    ############################################################################
#
#
#    ### Catch no given parameters ##############################################
#    if (!@ARGV && ${$paramsRef}{"NOPARAM"} && ref(${$paramsRef}{"NOPARAM"}) eq "CODE") {
#        print "\nWARNING: Parameters needed...\n\n";
#        &{${$paramsRef}{"NOPARAM"}}();
#    }
#    ############################################################################
#}



sub inflateBoxXY {
    $_[0]{'cooX'} *= $_[1];
    $_[0]{'cooY'} *= $_[1];
}



sub centerGroupXY {
    my $atomIdsRef = shift;
    my $coordsRef  = shift;
    my $boxRef     = shift;

    my %groupGeoCenter = getGeoCenterXY($atomIdsRef, $coordsRef);
    my %translVector = ('cooX' => ($$boxRef{'cooX'}*0.5-$groupGeoCenter{'cooX'}),
                        'cooY' => ($$boxRef{'cooY'}*0.5-$groupGeoCenter{'cooY'}));
    foreach (@{$atomIdsRef}) {
        $$coordsRef[$_]{'cooX'} += $translVector{'cooX'};
        $$coordsRef[$_]{'cooY'} += $translVector{'cooY'};
    }
}



sub getGeoCenterXY {
    my $atomIdsRef = shift;
    my $coordsRef  = shift;
    my %geoCenter  = ('cooX' => 0,
                      'cooY' => 0);

    foreach (@{$atomIdsRef}) {
        $geoCenter{'cooX'} += $$coordsRef[$_]{'cooX'};
        $geoCenter{'cooY'} += $$coordsRef[$_]{'cooY'};
    }
    $geoCenter{'cooX'} /= scalar(@$atomIdsRef);
    $geoCenter{'cooY'} /= scalar(@$atomIdsRef);
    return %geoCenter;
}



#sub getAtomIdList {
#    my $atomsRef = shift;
#    my $key      = shift;
#    my $needle   = shift;
#    my (@atomIdList, @restAtomIdList);
#
#    for (my $i=1; $i<@{$atomsRef}; $i++) {
#        $$atomsRef[$i]{$key} =~ /$needle/ ? push(@atomIdList, $i) : push(@restAtomIdList, $i);
#    }
#    return (\@atomIdList);
##    return (\@atomIdList, \@restAtomIdList);
#}



sub updateDelResnames {
    my $delResnamesRef  = shift;
    my $delAtomsDataRef = shift;
    my %tmpResid;

    for (my $i=0; $i<@{$delAtomsDataRef}; $i++) {
        next unless $$delAtomsDataRef[$i];
        next if $tmpResid{$$delAtomsDataRef[$i]{'uResId'}};
        $tmpResid{$$delAtomsDataRef[$i]{'uResId'}} = 1;
        $$delResnamesRef{$$delAtomsDataRef[$i]{'resName'}}++;
    }
}



sub renumAtoms {
    my $coordsRef = shift;
    my @newCoords;
    my $newAtomId = 1;
    my @renumMatrix;

    for (my $oldAtomId=0; $oldAtomId<@{$coordsRef}; $oldAtomId++) {
        next unless $$coordsRef[$oldAtomId]{'uResId'};
        $renumMatrix[$oldAtomId] = $newAtomId;
        $newCoords[$newAtomId++] = $$coordsRef[$oldAtomId];
    }
    return (\@newCoords, @renumMatrix);
}



#sub resname2moltype {
#    my $itpIncludeRef = shift;
#    my $resname       = shift;
#    my $topPath       = shift;
#    my $molType;
#
#    for (my $i=0; $i<@{$itpIncludeRef}; $i++) {
#        my $itpFile = -e ($topPath . $$itpIncludeRef[$i]) ? $topPath . $$itpIncludeRef[$i] : ($GMXLIB . "/" . $$itpIncludeRef[$i]);
#        print "Checking ITP file \"$itpFile\" for residue name $resname\n";
#        $molType = ITPFiles::resname2moltype($itpFile, $resname);
#        print "  Found residue name $resname in ITP file \"$itpFile\": moleculetype = \"$molType\"\n" if $molType;
#        return $molType if $molType;
#    }
#}



sub round {
    my ( $num, $prec ) = @_;
    return int( $num / $prec + 0.5 - ( $num < 0 ) ) * $prec;
}



sub getMax {
    return $_[0] if $_[0] > $_[1];
    return $_[1];
}



sub getMin {
    return $_[0] if $_[0] < $_[1];
    return $_[1];
}
################################################################################



################################################################################
### Lipid routines #############################################################
################################################################################
package Lipid;

sub analyzeBilayer {
    my $atomIdsRef   = shift;
    my $coordsRef    = shift;
    my $refAtomName  = shift;
    my @headAtomIds;
    my %lipidData    = ('zCenter' => 0);

    foreach (@{$atomIdsRef}) {
        next unless $$coordsRef[$_]{'atomName'};
        $lipidData{'zMin'} = $$coordsRef[$_]{'cooZ'} unless $lipidData{'zMin'};
        $lipidData{'zMax'} = $$coordsRef[$_]{'cooZ'} unless $lipidData{'zMax'};
        $lipidData{'zMin'} = $$coordsRef[$_]{'cooZ'} if $$coordsRef[$_]{'cooZ'} < $lipidData{'zMin'};
        $lipidData{'zMax'} = $$coordsRef[$_]{'cooZ'} if $$coordsRef[$_]{'cooZ'} > $lipidData{'zMax'};
        next unless $$coordsRef[$_]{'atomName'} =~ /$refAtomName/;
        $lipidData{'zCenter'} += $$coordsRef[$_]{'cooZ'};
        push(@headAtomIds, $_);
    }
    $lipidData{'zCenter'} /= scalar(@headAtomIds) if @headAtomIds;

    $lipidData{'leaflets'} = detectLeaflets(\@headAtomIds, $coordsRef, $lipidData{'zCenter'});

    return %lipidData;
}



sub detectLeaflets {
    my $headAtomIdsRef = shift;
    my $coordsRef      = shift;
    my $bilayerCenterZ = shift;
    my %leafletData;
    $leafletData{'upper'}{'nLipids'} = 0;
    $leafletData{'lower'}{'nLipids'} = 0;

    foreach (@{$headAtomIdsRef}) {
        if ($$coordsRef[$_]{'cooZ'} > $bilayerCenterZ) {
            $leafletData{'upper'}{'uResIds'}[ $$coordsRef[$_]{'uResId'} ] = 1;
            $leafletData{'upper'}{'headCooZAvg'} += $$coordsRef[$_]{'cooZ'};
            $leafletData{'upper'}{'nLipids'}++;
        }
        else {
            $leafletData{'lower'}{'uResIds'}[ $$coordsRef[$_]{'uResId'} ] = 1;
            $leafletData{'lower'}{'headCooZAvg'} += $$coordsRef[$_]{'cooZ'};
            $leafletData{'lower'}{'nLipids'}++;
        }
    }
    $leafletData{'upper'}{'headCooZAvg'} /= $leafletData{'upper'}{'nLipids'} if $leafletData{'upper'}{'nLipids'};
    $leafletData{'lower'}{'headCooZAvg'} /= $leafletData{'lower'}{'nLipids'} if $leafletData{'lower'}{'nLipids'};

    return (\%leafletData);
}
################################################################################



################################################################################
### Protein routines ###########################################################
################################################################################
#package Protein;
#
#sub analyzeProtein {
#    my $ndxDataRef      = shift;
#    my $protGroupId     = shift;
#    my $suGroupIdsRef   = shift;
#    my $coordsRef       = shift;
#    my $boxRef          = shift;
#    my $gridDelta       = shift;
#    my $zMin            = shift;
#    my $zMax            = shift;
#    my $extCavDetection = shift;
#    my $xyProfile       = shift;
#
#    my $gridRef;
#    my $gridsizeX = 0;
#    my $gridsizeY = 0;
#
#    my $gridAreaProtein;
#    my $gridAreaCavity;
#
#
#    ### Initialize grid ########################################################
#    ($gridRef, $gridsizeX, $gridsizeY) = iniGrid($boxRef, $gridDelta);
#    ############################################################################
#
#
#    ### Create protein grid ####################################################
#    $gridAreaProtein = protein2Grid($gridRef, $gridsizeX, $gridsizeY, $gridDelta, $$ndxDataRef[$protGroupId]{'atoms'}, $coordsRef, $zMin, $zMax);
#    ############################################################################
#
#
#    ### Connect subdomains to form a cavity ####################################
#    connectSubunits($gridRef, $gridsizeX, $gridsizeY, $gridDelta, $ndxDataRef, $suGroupIdsRef, $coordsRef, $zMin, $zMax) if $suGroupIdsRef;
#    ############################################################################
#
#
#    ### Detect the cavity ######################################################
#    $gridAreaCavity = detectCavity($gridRef, $gridsizeX, $gridsizeY, $gridDelta, , $extCavDetection);
#    ############################################################################
#
#
#    ### Write out the XY profile ###############################################
#    writeXyProfile ($xyProfile, $gridRef, $gridsizeX, $gridsizeY, $gridDelta) if $xyProfile;
#    ############################################################################
#
#    return ($gridAreaProtein, $gridAreaCavity, $gridRef);
#}
#
#
#
#sub iniGrid {
#    my $boxRef    = shift;
#    my $gridDelta = shift;
#
#    my @grid;
#    my $gridsizeX = int($$boxRef{'cooX'} / $gridDelta) + 1;
#    my $gridsizeY = int($$boxRef{'cooY'} / $gridDelta) + 1;
#
#    printf("      Create a %dx%d XY grid: 0%%\r", $gridsizeX, $gridsizeY) if $main::verbose;
#    for (my $x=0; $x<=$gridsizeX; $x++) {
#        printf("      Create a %dx%d XY grid: %d%%\r", $gridsizeX, $gridsizeY, 100*$x/$gridsizeX) if $main::verbose;
#        for (my $y=0; $y<=$gridsizeY; $y++) {
#            $grid[$x][$y] = 0;
#        }
#    }
#    printf("      Create a %dx%d XY grid: 100%%\n", $gridsizeX, $gridsizeY) if $main::verbose;
#    ############################################################################
#
#    return (\@grid, $gridsizeX, $gridsizeY);
#}
#
#
#
#sub protein2Grid {
#    my $gridRef     = shift;
#    my $gridsizeX   = shift;
#    my $gridsizeY   = shift;
#    my $gridDelta   = shift;
#    my $atomsIdsRef = shift;
#    my $coordsRef   = shift;
#    my $zMin        = shift;
#    my $zMax        = shift;
#
#    my $nAreas      = 0;
#    my $gridDelta2  = $gridDelta*$gridDelta;
#
#    print "      Mapping atoms to the grid: 0.000%\r" if $main::verbose;
#    foreach (@{$atomsIdsRef}) {
#        next unless $$coordsRef[$_]{'cooZ'};
#        next if $$coordsRef[$_]{'cooZ'} > $zMax;
#        next if $$coordsRef[$_]{'cooZ'} < $zMin;
#        printf("      Mapping atoms to the grid: %d%% (protein area %.4f nm^2)\r", (++$$coordsRef[$_]{'atomNum'})/scalar(@{$atomsIdsRef})*100, $nAreas*$gridDelta*$gridDelta) if $main::verbose;
#
#        my $element = substr($$coordsRef[$_]{'atomName'}, 0, 1);
#        my $radius  = PTE::getRadius($element);
#        my $radius2 = $radius * $radius;
#
#        my $tmpGridX = int($$coordsRef[$_]{'cooX'} / $gridDelta);
#        my $tmpGridY = int($$coordsRef[$_]{'cooY'} / $gridDelta);
#        my $subrange = int($radius / $gridDelta);
#
#        for (my $x=getMax($tmpGridX-$subrange, 0); $x<=getMin($tmpGridX+$subrange, $gridsizeX-1); $x++) {
#            for (my $y=getMax($tmpGridY-$subrange, 0); $y<=getMin($tmpGridY+$subrange, $gridsizeY-1); $y++) {
#                my $dx = $$coordsRef[$_]{'cooX'} - $x * $gridDelta;
#                my $dy = $$coordsRef[$_]{'cooY'} - $y * $gridDelta;
#                my $dist2 = $dx*$dx + $dy*$dy;
#                next if $dist2 > $radius2;
#                ++$nAreas unless $$gridRef[$x][$y];
#                $$gridRef[$x][$y] = 1;
#            }
#        }
#    }
#    printf("      Mapping atoms to the grid: 100%% (protein area %.4f nm^2)\n", $nAreas*$gridDelta2) if $main::verbose;
#
#    return ($nAreas*$gridDelta2);
#}
#
#
#
#sub connectSubunits {
#    my $gridRef       = shift;
#    my $gridsizeX     = shift;
#    my $gridsizeY     = shift;
#    my $gridDelta     = shift;
#    my $ndxDataRef    = shift;
#    my $suGroupIdsRef = shift;
#    my $coordsRef     = shift;
#    my $zMin          = shift;
#    my $zMax          = shift;
#
#    my @gridSuGeoCenter;
#
#
#    print "      Detecting subunit connections: 0.0000 nm^2 possible\r" if $main::verbose;
#
#    ### Calculate center of mass of each protein subunit on the grid ###########
#    foreach my $suGroupId (@{$suGroupIdsRef}) {
#        my %tmpSum = ('cooX' => 0, 'cooY' => 0);
#        my $nAtoms = 0;
#        foreach (@{$ndxData[$suGroupId]{'atoms'}}) {
#            next unless $$coordsRef[$_]{'cooZ'};
#            next if $$coordsRef[$_]{'cooZ'} > $zMax;
#            next if $$coordsRef[$_]{'cooZ'} < $zMin;
#            $tmpSum{'cooX'} += $$coordsRef[$_]{'cooX'};
#            $tmpSum{'cooY'} += $$coordsRef[$_]{'cooY'};
#            $nAtoms++;
#        }
#        next unless $nAtoms;
#        $tmpSum{'cooX'} = int(($tmpSum{'cooX'}/$nAtoms)/$gridDelta);
#        $tmpSum{'cooY'} = int(($tmpSum{'cooY'}/$nAtoms)/$gridDelta);
#        push(@gridSuGeoCenter, \%tmpSum);
#    }
#    ############################################################################
#
#
#    ### Print the geometrical center of each subunit ###########################
#    for (my $i=0; $i<@gridSuGeoCenter; $i++) {
#        getGridSlope($gridRef, $gridSuGeoCenter[$i-1], $gridSuGeoCenter[$i]);
#    }
#    ############################################################################
#}
#
#
#
## This does not work for all cases. Replace this routine by a vector based method.
#sub getGridSlope {
#    my $gridRef   = shift;
#    my $vecARef   = shift;
#    my $vecBRef   = shift;
#
#    my $slope = ($$vecARef{'cooY'} - $$vecBRef{'cooY'}) / ($$vecARef{'cooX'} - $$vecBRef{'cooX'});
#    my $xMin;
#    my $xMax;
#    my $yMin;
#
#    if ($$vecARef{'cooX'} < $$vecBRef{'cooX'}) {
#        $xMin = $$vecARef{'cooX'};
#        $yMin = $$vecARef{'cooY'};
#
#        $xMax = $$vecBRef{'cooX'};
#    }
#    else {
#        $xMin = $$vecBRef{'cooX'};
#        $yMin = $$vecBRef{'cooY'};
#
#        $xMax = $$vecARef{'cooX'};
#    }
#
#    if ($slope > 3) {
#        $xMin = $$vecARef{'cooY'};
#        $yMin = $$vecARef{'cooX'};
#
#        $xMax = $$vecBRef{'cooY'};
#    }
#
#    #if ($$vecARef{'cooX'} > $$vecBRef{'cooX'}) {
#    #    $xMax = $$vecARef{'cooX'};
#    #}
#    #else {
#    #    $xMax = $$vecBRef{'cooX'};
#    #}
#
#    for (my $x=$xMin; $x<=$xMax; $x++) {
#        $yMin += $slope;
#        my $intY = main::round($yMin, 1);
#        next if $$gridRef[$x][$intY];
#        $$gridRef[$x][$intY] = 3;
#
#        $$gridRef[$x+1][$intY] = 3;
#        $$gridRef[$x-1][$intY] = 3;
#        $$gridRef[$x][$intY+1] = 3;
#        $$gridRef[$x][$intY-1] = 3;
#        $$gridRef[$x+1][$intY+1] = 3;
#        $$gridRef[$x+1][$intY-1] = 3;
#        $$gridRef[$x-1][$intY+1] = 3;
#        $$gridRef[$x-1][$intY-1] = 3;
#    }
#}
#
#
#
#sub detectCavity {
#    my $gridRef         = shift;
#    my $gridsizeX       = shift;
#    my $gridsizeY       = shift;
#    my $gridDelta       = shift;
#    my $extCavDetection = shift;
#
#    my $nCavityAreas  = 0;
#    my $gridDelta2    = $gridDelta*$gridDelta;
#
#
#    print "      Detecting protein internal cavities: 0.0000 nm^2 possible\r" if $main::verbose;
#    ### Invert protein grid (set all empty grid points to 2 (poss. cav.)) ######
#    for (my $x=0; $x<=$gridsizeX; $x++) {
#        for (my $y=0; $y<=$gridsizeY; $y++) {
#            next if $$gridRef[$x][$y]; # Next loop if this grid point is defined as protein.
#            $$gridRef[$x][$y] = 2;
#            $nCavityAreas++;
#        }
#        printf("      Detecting protein internal cavities: %.4f nm^2 possible      \r", $nCavityAreas*$gridDelta2) if $main::verbose;
#    }
#    ############################################################################
#
#
#    ### Washing out the cavities ###############################################
#    for (my $i=0; $i<2; $i++) {
#        my $foundExcl1 = 1;
#        my $foundExcl2 = 1;
#        for (my $x=0; $x<=$gridsizeX; $x++) {
#            my $neighbCellX = $x - 1;
#            $nCavityAreas -= $foundExcl1 = washingOutY($x, $neighbCellX, $gridsizeY, $gridRef, $extCavDetection);
#            printf("      Detecting protein internal cavities: %.4f nm^2 possible      \r", $nCavityAreas*$gridDelta2) if $main::verbose;
#        }
#
#        for (my $x=$gridsizeX; $x>=0; $x--) {
#            my $neighbCellX = $x + 1;
#            $nCavityAreas -= $foundExcl2 = washingOutY($x, $neighbCellX, $gridsizeY, $gridRef, $extCavDetection);
#            printf("      Detecting protein internal cavities: %.4f nm^2 possible      \r", $nCavityAreas*$gridDelta2) if $main::verbose;
#        }
#        $i = 0 if $foundExcl1 || $foundExcl2;
#    }
#    ############################################################################
#    printf("      Detecting protein internal cavities: %.4f nm^2 detected      \n", $nCavityAreas*$gridDelta2) if $main::verbose;
#
#    return ($nCavityAreas*$gridDelta2);
#}
#
#
#
#sub washingOutY {
#    my $x               = shift;
#    my $neighbCellX     = shift;
#    my $gridsizeY       = shift;
#    my $gridRef         = shift;
#    my $extCavDetection = shift;
#
#    my $delAreas        = 0;
#
#    for (my $y=0; $y<=$gridsizeY; $y++) { # SW -> NE.
#        my $neighbCellY = $y - 1;
#        next unless $$gridRef[$x][$y] == 2; # Next if this grid point is not a possible cavity.
#        unless ($$gridRef[$x][$neighbCellY]) { # If the neighbored cell is not protein or not a possible cavity...
#            $$gridRef[$x][$y] = 0; # ...exclude this from a possible cavity.
#            $delAreas++;
#            next;
#        }
#        unless ($$gridRef[$neighbCellX][$y]) {
#            $$gridRef[$x][$y] = 0;
#            $delAreas++;
#            next;
#        }
#
#        next unless $extCavDetection;
#        unless ($$gridRef[$neighbCellX][$neighbCellY]) {
#            $$gridRef[$x][$y] = 0;
#            $delAreas++;
#        }
#    }
#
#    for (my $y=$gridsizeY; $y>=0; $y--) { # NW -> SE.
#        my $neighbCellY = $y + 1;
#        next unless $$gridRef[$x][$y] == 2;
#        unless ($$gridRef[$x][$neighbCellY]) {
#            $$gridRef[$x][$y] = 0;
#            $delAreas++;
#            next;
#        }
#        unless ($$gridRef[$neighbCellX][$y]) {
#            $$gridRef[$x][$y] = 0;
#            $delAreas++;
#            next;
#        }
#
#        next unless $extCavDetection;
#        unless ($$gridRef[$neighbCellX][$neighbCellY]) {
#            $$gridRef[$x][$y] = 0;
#            $delAreas++;
#        }
#    }
#
#    return $delAreas;
#}
#
#
#
#sub writeXyProfile {
#    my $xyProfile = shift;
#    my $gridRef   = shift;
#    my $gridsizeX = shift;
#    my $gridsizeY = shift;
#    my $gridDelta = shift;
#
#    printf("      Writing out protein xy grid profile: 0%%\r") if $main::verbose;
#    open(XYPROFILE, ">$xyProfile") || die "ERROR: Cannot open profile output file \"$xyProfile\": $!\n";
#    for (my $x=0; $x<=$gridsizeX; $x++) {
#        my $tmpX = $x*$gridDelta; # Performance.
#        for (my $y=0; $y<=$gridsizeY; $y++) {
#            printf(XYPROFILE "%f %f %d\n", $tmpX, $y*$gridDelta, $$gridRef[$x][$y]) if $$gridRef[$x][$y];
#        }
#        printf("      Writing out protein xy grid profile: %d%%\r", 100*$x/$gridsizeX) if $main::verbose;
#    }
#    close XYPROFILE;
#    printf("      Writing out protein xy grid profile: 100%%\n") if $main::verbose;
#}
#
#
#
#sub getMax {
#    return $_[0] if $_[0] > $_[1];
#    return $_[1];
#}
#
#
#
#sub getMin {
#    return $_[0] if $_[0] < $_[1];
#    return $_[1];
#}
################################################################################



################################################################################
### Periodic Table of the Elements routines ####################################
################################################################################
#package PTE;
#
#sub getRadius {
#    my %atomRadius = ('H'  => 0.110,
#                      'He' => 0.140,
#                      'Li' => 0.182,
#                      'Be' => 0.153,
#                      'B'  => 0.192,
#                      'C'  => 0.170,
#                      'N'  => 0.155,
#                      'O'  => 0.152,
#                      'F'  => 0.147,
#                      'Ne' => 0.154,
#                      'Na' => 0.227,
#                      'Mg' => 0.173,
#                      'Al' => 0.184,
#                      'Si' => 0.210,
#                      'P'  => 0.180,
#                      'S'  => 0.180,
#                      'Cl' => 0.175,
#                      'Ar' => 0.188,
#                      'K'  => 0.275,
#                      'Ca' => 0.231,
#                      'Ni' => 0.163,
#                      'Cu' => 0.140,
#                      'Zn' => 0.139,
#                      'Ga' => 0.187,
#                      'Ge' => 0.211,
#                      'As' => 0.185,
#                      'Se' => 0.190,
#                      'Br' => 0.185,
#                      'Kr' => 0.202,
#                      'Rb' => 0.303,
#                      'Sr' => 0.249,
#                      'Pd' => 0.163,
#                      'Ag' => 0.172,
#                      'Cd' => 0.158,
#                      'In' => 0.193,
#                      'Sn' => 0.217,
#                      'Sb' => 0.206,
#                      'Te' => 0.206,
#                      'I'  => 0.198,
#                      'Xe' => 0.216,
#                      'Cs' => 0.343,
#                      'Ba' => 0.268,
#                      'Pt' => 0.175,
#                      'Au' => 0.166,
#                      'Hg' => 0.155,
#                      'Tl' => 0.196,
#                      'Pb' => 0.202,
#                      'Bi' => 0.207,
#                      'Po' => 0.197,
#                      'At' => 0.202,
#                      'Rn' => 0.220,
#                      'Fr' => 0.348,
#                      'Ra' => 0.283,
#                      'U'  => 0.186);
#
#    return $atomRadius{$_[0]} if $atomRadius{$_[0]};
#    return 0;
#
## Atomic van der Waals radii in nm.
##    1) R. Scott Rowland, Robin Taylor: Intermolecular Nonbonded Contact
##       Distances in Organic Crystal Structures: Comparison with Distances
##       Expected from van der Waals Radii. In: J. Phys. Chem. 1996, 100,
##       S. 7384–7391, doi:10.1021/jp953141+.
##    2) A. Bondi: van der Waals Volumes and Radii. In: J. Phys. Chem. 1964, 68,
##       S. 441-451, doi:10.1021/j100785a001.
##    3) Manjeera Mantina, Adam C. Chamberlin, Rosendo Valero, Christopher J.
##       Cramer, Donald G. Truhlar: Consistent van der Waals Radii for the Whole
##       Main Group. In: J. Phys. Chem. A. 2009, 113, S. 5806–5812,
##       doi:10.1021/jp8111556.
#}
################################################################################




################################################################################
### PDBFiles  ##################################################################
################################################################################
package PDBFiles;

sub readPdb {
    my $pdbFile = shift;
    my %pdbData;

    print "  ---------------------------------\n  Read PDB file \"$pdbFile\"...\r";
    open(PDBFILE, "<$pdbFile") || die "ERROR: Cannot open PDB file \"$pdbFile\": $!\n";
    readCoords(\*PDBFILE, \%pdbData);
    close(PDBFILE);
    print "  Read PDB file \"$pdbFile\": Finished\n  ---------------------------------\n\n";

    return %pdbData;
}



sub readCoords {
    my $fileHandle = shift;
    my $pdbDataRef = shift;
    my $atomId     = 0;
    my $uResId     = 0; # The unique residue ID.

    while (<$fileHandle>) {
        chomp($_);
        $$pdbDataRef{'atoms'}[++$atomId] = getAtomdata($_) if ($_ =~ /^ATOM\s+/);
        $$pdbDataRef{'atoms'}[$atomId]{'uResId'} = $uResId++ unless ($_ =~ /^ATOM\s+/);
        print "    Read atom data:  $atomId\r" if $main::verbose;
    }
    return 1;
}



sub getAtomdata {
    my $atomStr = shift;
    my $strLen  = length($atomStr);
    my %atomData;

    $atomData{'atomNum'}    = checkSubstr($atomStr, $strLen, 6, 5);
    $atomData{'atomName'}   = checkSubstr($atomStr, $strLen, 12, 4);
    $atomData{'altLoc'}     = checkSubstr($atomStr, $strLen, 16, 1);
    $atomData{'resName'}    = checkSubstr($atomStr, $strLen, 17, 3);
    $atomData{'chainID'}    = checkSubstr($atomStr, $strLen, 21, 1);
    $atomData{'resId'}      = checkSubstr($atomStr, $strLen, 22, 4);
    $atomData{'iCode'}      = checkSubstr($atomStr, $strLen, 26, 1);
    $atomData{'cooX'}       = checkSubstr($atomStr, $strLen, 30, 8);
    $atomData{'cooY'}       = checkSubstr($atomStr, $strLen, 38, 8);
    $atomData{'cooZ'}       = checkSubstr($atomStr, $strLen, 46, 8);
    $atomData{'occupancy'}  = checkSubstr($atomStr, $strLen, 54, 6);
    $atomData{'tempFactor'} = checkSubstr($atomStr, $strLen, 60, 6);
    $atomData{'element'}    = checkSubstr($atomStr, $strLen, 76, 2);
    $atomData{'charge'}     = checkSubstr($atomStr, $strLen, 78, 2);

    return \%atomData;
}



sub checkSubstr {
    my $str       = shift;
    my $strLen    = shift;
    my $start     = shift;
    my $substrLen = shift;
    my $substr    = '';

    if ($strLen >= ($start+$substrLen)) {
        $substr = substr($str, $start, $substrLen);
        $substr =~ s/\s//g;
    }
    return $substr;
}



sub writePdb {
    my $pdbFile    = shift;
    my $pdbDataRef = shift;

    $pdbFile .= ".pdb" unless $pdbFile =~ /\.pdb$/;

    open(PDBFILE, ">$pdbFile") || die "ERROR: Cannot open output PDB file ($pdbFile): $!\n";
    writeCoords(\*PDBFILE, $pdbDataRef);
    close(PDBFILE);
}



sub writeCoords {
    my $fileHandle = shift;
    my $pdbDataRef = shift;

    my $atomId     = 1;

    my %defaultAtom = ("atomNum"    => 0,
                       "atomName"   => "X",
                       "altLoc"     => "",
                       "resName"    => "RES",
                       "chainId"    => 0,
                       "resId"      => 0,
                       "iCode"      => "",
                       "cooX"       => 0.0,
                       "cooY"       => 0.0,
                       "cooZ"       => 0.0,
                       "occupancy"  => 0.0,
                       "tempFactor" => 0.0,
                       "element"    => "X",
                       "charge"     => "0");

    foreach (@{$$pdbDataRef{'atoms'}}) {
        next unless $$_{'uResId'};

        ### Fill empty fields with standard parameters #########################
        foreach my $key (keys %defaultAtom) {
            $$_{$key} = $defaultAtom{$key} unless defined($$_{$key});
        }
        ########################################################################
        $$_{'element'} = substr($$_{'atomName'}, 0, 1) if $$_{'element'} eq 'X';

        printf($fileHandle "ATOM  %5d %4s%1s%4s%1s%4d%1s   %8.3f%8.3f%8.3f%6.2f%6s          %2s%2s\n",
            (($atomId++)%100000), $$_{'atomName'}, $$_{'altLoc'}, $$_{'resName'}, $$_{'chainId'}, ($$_{'resId'}%100000),  $$_{'iCode'}, ($$_{'cooX'}*10), ($$_{'cooY'}*10), ($$_{'cooZ'}*10), $$_{'occupancy'}, substr(sprintf("%6.2f", $$_{'tempFactor'}), 0, 6), $$_{'element'}, $$_{'charge'});
    }
}
################################################################################



################################################################################
### TOPFiles  ##################################################################
################################################################################
package TOPFiles;

sub readTop {
    my $topFile = shift;
    my %topData;

    print "  ---------------------------------\n  Read TOP file \"$topFile\"...\r";
    open(TOPFILE, "<$topFile") || die "ERROR: Cannot open TOP file \"$topFile\": $!\n";
    readData(\*TOPFILE, \%topData);
    close(TOPFILE);
    print "  Read TOP file \"$topFile\": Finished\n  ---------------------------------\n\n";

    return %topData;
}



sub readData {
    my $fileHandle = shift;
    my $topDataRef = shift;
    my $molSwitch  = 0;

    while (<$fileHandle>) {
        chomp($_);
        push(@{$$topDataRef{'include'}}, $1) if $_ =~ /^\s*#include\s+"(.+)"/;
        $molSwitch = 0 if $molSwitch && $_ =~ /^\s*\[.+\]"/;
        $$topDataRef{'molecules'}{$1} = $2 if $molSwitch && $_ =~ /^\s*(\w+)\s+(\d+)/;
        $molSwitch = 1 if $_ =~ /^\s*\[ molecules \]/;
    }
    printf("\n    Read files for inclusion: %d\n", scalar @{$$topDataRef{'include'}}) if $main::verbose;
    printf("    Number of molecule types: %d\n", scalar keys %{$$topDataRef{'molecules'}}) if $main::verbose;
}



sub updateTop {
    my $topFile = shift;
    my $molType = shift;
    my $molNum  = shift;

    my $backupFile = main::fileBackup($topFile);
    open(BUTOPFILE, "<$backupFile") || die "ERROR: Cannot open TOP file \"$topFile\": $!\n";
    open(TOPFILE, ">$topFile") || die "ERROR: Cannot open TOP file \"$topFile\": $!\n";
    while (<BUTOPFILE>) {
        $_ = sprintf("%-14s %6d\n", $molType, $molNum) if $_ =~ /^\s*$molType/;
        print TOPFILE $_;
    }
    close(TOPFILE);
    close(BUTOPFILE);
}
################################################################################


###############################################################################
### APLFiles  (APL = area per lipid) ##########################################
###############################################################################
package APLFiles;

sub writeApl {
    my $aplFile    = shift;
    my $aplDataRef = shift;

    my $isNewFile = (-e $aplFile) ? 0 : 1;
    open(APLFILE, ">>$aplFile") || die "ERROR: Cannot open area per lipid output file \"$aplFile\": $!\n";
    print APLFILE "# Group-ID | Group-Name | Area/lipid (upper leaflet) | Area/lipid (lower leaflet) | Apl(average [both leaflets])\n" if $isNewFile;
    printf APLFILE ("%10d | %10s", $$aplDataRef{'groupId'}, $$aplDataRef{'groupName'});
    printf APLFILE (" | %f", $$aplDataRef{'upper'}) if $$aplDataRef{'upper'};
    printf APLFILE (" | %f", $$aplDataRef{'lower'}) if $$aplDataRef{'lower'};
    printf APLFILE (" | %f", $$aplDataRef{'average'}) if $$aplDataRef{'average'};
    print APLFILE "\n";
    close(APLFILE);
}
################################################################################
